/***********************************************************************
Compile       : gcc -O3 -o matrix matrix.c ../design/design.c ../design/utils.c dynamics.c -I../design
Usage         : ./matrix mode mode1 mode2 irflag N L inputfile 
                (./matrix 0 0 0 0 11 12 Budding.txt)
mode          : 0 normal ; 1 cycled ; 2 triggled cycled
mode1         : 0 disallow self-activation ; 1 allow self_activation
mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II
irflag        : 1 irreducible solution only
N             : number of the nodes
L             : number of the steps        
Output        : designability matrix
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include "design.h"
#include "dynamics.h"

void print_help ()
{
	printf ("Usage         : ./matrix mode mode1 mode2 irflag N L inputfile\n");
	printf ("                (./matrix 0 0 0 0 11 12 Budding.txt)\n");
	printf ("mode          : 0 normal ; 1 cycled ; 2 triggled cycled\n");
	printf ("mode1         : 0 disallow self-activation ; 1 allow self_activation\n");
	printf ("mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II\n");
	printf ("irflag        : 1 irreducible solution only\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n\n");
	printf ("Output        : designability matrix\n");
}

void init_design (struct designability *design, int argc, char *argv[])
{
	int N, L, i;

	design->mode = atoi (argv[1]);
	design->mode1 = atoi (argv[2]);
	design->mode2 = atoi (argv[3]);
	design->irflag = atoi (argv[4]);
	design->bdflag = 0;
	design->dgflag = 0;
	design->pnflag = 0;
	design->svflag = 1;
	design->N = atoi (argv[5]);
	design->L = atoi (argv[6]);

	design->N += design->mode1;
	N = design->N;
	L = design->L;

	/* the degree of irreducible solution is zero */
	if (design->irflag)
		design->dgflag = 0;

	design->pathway = calloc (sizeof (int), N * L);
	design->NS2 = calloc (sizeof (int), N);
	design->min_link = calloc (sizeof (int), N);
	design->links_dist = calloc (sizeof (int), N * (N + 1));
	design->Links_dist = calloc (sizeof (double), N * N + 1);
	design->degree_dist = calloc (sizeof (int), N * (2 * N + 1));
	design->Degree_dist = calloc (sizeof (double), 2 * N * N + 1);
	design->solution = calloc (sizeof (struct solution), N);

	for (i = 0; i < N; i++)
	{
		init_solution (&(design->solution[i]));
	}
}

void free_design (struct designability *design)
{
	int i;

	free (design->pathway);
	free (design->NS2);
	free (design->min_link);
	free (design->links_dist);
	free (design->Links_dist);
	free (design->degree_dist);
	free (design->Degree_dist);
	free (design->solution);

	for (i = 0; i < design->N; i++)
	{
		free_solution (&(design->solution[i]));
	}
}

void print_solution (struct designability *design)
{
	int i, j, k, N;

	N = design->N;
	if (design->svflag)
	{
		for (i = design->mode1; i < N; i++)
		{
			printf ("Total Number of Soultion for Node %2d : %d\n", i - design->mode1, design->NS2[i]);
			for (j = 0; j < design->solution[i].curr; j++)
			{
				printf ("Node %2d :", i - design->mode1);
				for (k = design->mode1; k < N; k++)
				{
					if (design->mode1 == 1 && k == i && design->solution[i].edge[j * N] == 1)
						printf (" g");
					else if (design->solution[i].edge[j * N + k] == 0)
						printf (" n");
					else if (design->solution[i].edge[j * N + k] == 1)
						printf (" g");
					else if (design->solution[i].edge[j * N + k] == -1)
						printf (" r");
				}
				printf ("\n");
			}
		}
	}
}

int main (int argc, char *argv[])
{
	struct designability design;

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	init_design (&design, argc, argv);

	FILE *fp;
	char buf[256];
	int i, j, k, N, L;

	N = design.N;
	L = design.L;

	fp = fopen (argv[7], "r");

	fgets (buf, 256, fp);
	design.activity = 0;
	for (i = 0; i < L; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (design.mode1 == 1 && j == 0)
			{
				design.pathway[i * N + j] = 1;
				continue;
			}
			fscanf (fp, "%d", &(design.pathway[i * N + j]));
			if (i <= L - 3 || (i == L - 2 && design.mode != 2))
				design.activity += design.pathway[i * N + j];
		}
		fgets (buf, 256, fp);
	}

	if (design.mode == 2)
		design.activity --;

	/* Get the number of solution using enumerating algorithm */
	calc_design (&design);

	fclose (fp);

	int M;
	M = 1;
	for (i = design.mode1; i < N; i++)
	{
		M *= 2;
	}

	int state[N], state2[N], count[N];
	double NS[M][M];

	for (i = 0; i < M; i++)
	{
		unpack (i, N, state);
		if (design.mode1 == 1)
			state[0] = 1;
		for (j = design.mode1; j < N; j++)
		{
			count[j] = 0;
			for (k = 0; k < design.solution[j].curr; k++)
			{
				 count[j] += next_step (N, state, j, design.solution[j].edge + k * N);
			}
		}

		for (j = 0; j < M; j++)
		{
			unpack (j, N, state2);
			NS[i][j] = 1.0;
			for (k = design.mode1; k < N; k++)
			{
				if (state2[k] == 1)
					NS[i][j] *= count[k];
				else
					NS[i][j] *= (design.solution[k].curr - count[k]);
			}
			printf ("%04d -> %04d %lg\n", i, j, NS[i][j] / design.NS);
		}
	}

	fp = fopen ("DATA.DAT", "w");
	fwrite (NS, sizeof (double), M * M, fp);
	fclose (fp);

	free_design (&design);

	exit (0);
}
